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1. Introduction. 

Since the pioneering work of Poincare |TJ in celestial mechanics wherein 
the mathematical basis of deterministic chaos was established, the study of 
homoclinic phenomena has been the focus of increasing attention not only 
in celestial mechanics but in many different branches of physics, chemistry 
and biology [[J. Due to its universality and intrinsic mathematical interest, 
models in which unstable periodic orbits (UPOs) are subjected to small pe- 
riodic perturbations has been one of the main paradigms of deterministic 
chaos ||. An analytical tool to study these models is the Melnikov function 
that describes the transversal distance between the unstable and stable man- 
ifolds emanating from an UPO. Its isolated odd zeros indicate the crossing 
of these manifolds, hence the onset of chaos The use of Melnikov method 
in connection to Kolmogorov-Arnold-Moser (KAM) theory for a large class 
of perturbations was studied by Holmes and Marsden H, see also M. 

Examples of applications of the Melnikov method in gravitation are: 
the motion of a nearly symmetric heavy top J/J, the motion of particles 
in a two-dimensional Stackel potential separable with cosm<^ perturbations 
(m = 1, 2, ...) and other perturbations [Q, extensions of the former work for 
perturbation of the three-dimensional Stackel potentials H] the gravitational 



collapse of cosmological models [ItJ, the study of orbits around a black hole 
perturbed by either gravitational radiation O, |T2[ or an external quadrupo- 
lar shell JLJ 

In this paper we consider the equatorial motion of a particle moving in a 
potential modeled by an inverse square law plus a quadrupole-like term. This 
potential describes more realistically the gravitational force of attraction of 
stars (like the sun) on planets and galaxy bulges that to a certain degree of 
approximation can be represented by a monopole plus a quadrupole on stars 
outside the galaxy. Also, this potential arises naturally in General Relativity 
in the study of test particles in a spherically symmetric attraction center 
(Schwarzschild black hole), the "quadrupole" being associated to the test 
particle angular momentum, see for instance |14 . 
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Next, we consider a class of periodic perturbations of this integrable 
planar motion. The perturbations can be the result of periodic motion of 
mass distributions either inside or outside the orbit of interest. The fact that 
we were able to integrate in a closed form the equation of motion of the UPO 
is used to reduce the analysis of the Melnikov function to the study of simple 
graphics. We find that the perturbations will originate chaos in a variety of 
situations and range of parameters. Finally, we argue that chaotic behavior 
is generic for such a perturbations. 

2. The homoclinic orbit. 

We shall consider the orbit of a particle of mass m moving on a plane 
under the influence of a force described by a potential with a monopolar term 
plus a quadrupole-like contribution not necessarily small, 

V = -km/R - Qm/R 3 , (1) 

where k = GM, M is the mass of the attraction center, and Q is the 
quadrupolar strength. In astronomy, Q = GMJ 2 Rl/2, with R Q being a linear 
measure of the size of the central body and J 2 the distortion parameter, see 
for instance Hl4| . All the undefined symbols have their usual meaning. Note 



that for a body with reflection symmetry the octupolar term is null (J3 = 0). 

In General Relativity, Q is proportional to the square of the angular 
momentum of the test particle [|14[ . 



The Lagrangian of the particle is 



2 



Kit) +R UJ + R + 



(2) 



Since is a cyclic variable we have the conservation law h = R 2 d<j)/dt. Hence, 
the motion of the particle is determined by the effective Hamiltonian 

H = ^- + V eff , (P = mdR/dt), (3) 
V ^ = m {2R*-R-W- (4) 



3 



The stationary values of V e ff are R = ± ^1 — 12/3) with /3 = kQ/h 4 . 
The minus (plus) sign correspond to the relative maximum (minimum). Note 
that (3 is a non dimensional parameter. We shall denote by R un the UPO 
position given by the stationary value with minus sign (the maximum); we 
find that V e ff(R un ) = + 2y/l — (3). Since, for large R, the potential 

V e ff ~ — -| < 0, we need V e ff(R un ) < to have a periodic motion of the 
particle limited by R un and R m , the last being the turning point. Hence the 
parameter (3 is limited by 1/16 = 0.0625 < (3 < 1/12 ~ 0.0833. 

It is convenient to work with dimensionless quantities. We shall redefine 
our variables using the constants of the problem in the following way: r = 
kR/h 2 , t = k 2 t/h 3 } e = h 2 E/mk 2 , where E = H is the total energy of the 
particle. We recall that these quantities are proportional, respectively, to 
the radius, the period, and the energy of a particle of mass m moving in a 
circular orbit with areal velocity h around an attraction center characterized 
by k. The energy integral can be written as 



where the overdot indicates derivation with respect to the time parameter r. 
In these new units, R un is written as r un = (1 — y/1 — 12/3)/2. For f3 in the 
interval 1/16 < (3 < 1/12 we have 1/4 < r un < 1/2. In Fig. 1 we present 
a graphic of the potential @ for the values of (3 = 0.068 (top curve), 0.072, 
and 0.078 (bottom curve), the maximum value is located at r un . A particle 
with energy e = f e //(r un ) will describe either the UPO itself (r = r un , r = 0) 
or the homoclinic orbit tending to the UPO as r — > ±oo. This homoclinic 
orbit is enclosed by circles of radii r un and r m = 2r un (l — r un )/(Ar un — 1). 
Now, by making e = v e ff(r un ) we can write (|5|) in the equivalent form 



2e = r 2 + 2v eff , 

2v eff = 1/r 2 - 2/r - 2/3/r 3 , 



(5) 
(6) 



f 



(7) 



(r 



±W/3, 



with 




(8) 
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The differential equation (|7|) admits the quadrature 



T WpT = Jr(r m - r) + (r m + 2r un ) arctan 



ir m — r 



2r 3/2 

H — m ar ct anh , 



y/fm Tun \ ('"m Tun)f 



[r m - r r„ 



(9) 



We have chosen the constant of integration to have r = at r = r m . A 
graphic r(r) of the positive branch of (^|) is pictured in Fig. 2 for values of 
the parameter (3 = 0.064 (top curve), 0.065, 0.066, 0.067, and 0.068 (bottom 
curve). We see that the particle takes a finite time to travel from r m to the 
vicinity of r un and then an infinite time to arrive (depart) to (from) r un itself, 
wherein is located the UPO. 

3. The Melnikov method. 

2 

Let us consider an integrable Hamiltonian H = ^—+V(q), which admits 
at least one UPO with the corresponding homoclinic orbit, and also a small 
periodic perturbation described by the Hamiltonian function rjHi(q,p,t). 
Then the transverse distance in phase space between the unstable and the 
stable manifolds emanating from the UPO is proportional to 

/oo 
{H , H x }dt, (10) 
-oo 

where the integral is taken along the unperturbed homoclinic orbit, to is 
an arbitrary initial time, and {/, g} are the usual Poisson brackets, see for 
instance 0. If there is an intersection for some to, i.e., an isolated odd zero 
of M(to), then there will be one for every t . This infinitely many crossings 
of manifolds will produce a tangle that is the signature of homoclinic chaos 

B|]- 

The unperturbed system is Eqs. (§) and (|]) with H = e, and (q,p) = 
(r, r). For our purposes, it is enough to consider periodic perturbations of 
the particular form H\ = f(r) cos(f2r). Perturbations of this form are repre- 
sentative of a very general situation: take the potential of a mass distribution 
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moving periodically and make a Fourier expansion in time and a multipo- 
lar expansion in the space variables. Next, consider approximate axial and 
reflection symmetry. Therefore, the generic term of the series expansion in 
the plane of the UPO (6 = 7r/2) will look like Hi above with /(r) a positive 
(negative) powers of r for mass distributions outside (inside) the homoclinic 
orbit. We find for the Melnikov function (|T0|). 

In our case, for the homoclinic orbit (P), we get 

M(t ) = -2K{ti) sin(fir ), (12) 

K{Q) = f m sin[Qr(r)]dr, (13) 

J r un dr 

where in the integrand r(r) means the positive branch of (||). Note that in 
(|12"D does not appear a term proportional to cos(f2r ), its coefficient being 
null due to the fact that the homoclinic orbit has reflection symmetry with 
respect to the r-axis, and cos(f2r) is an even function in the r variable. Thus 
the Melnikov function will have the required zeros as long as K(Q) ^ 0. The 
coordinate change r — > r, instead of the usual r — ► 0, allows us to pass from 
an infinite interval to a finite one in (|HJ). This observation will play a key 
role in the evaluation of the function K(Q). 

4. Particular cases. 

Firstly, we shall consider the function K(Q) for exterior perturbations 
of the type f(r) = r n with n = 1,2, .... For a particle moving around the 
attraction center, they model periodic perturbations due to mass distribu- 
tions beyond the orbit of the particle, n = 1 represents dipolar contributions, 
n = 2 quadrupolar, etc. Static quadrupolar and octopolar exterior contribu- 
tions, inspired in the well known Henon-Heiles model, give rise also to chaotic 
behavior in General Relativity 



Now, we shall study the integrand of K(Q) formed by the product of an 
oscillating function by a polynomial. Near r = r un this oscillation is a rapid 
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one. To better understand this behavior we display in Fig. 3 a graphic of 
sin[fir(r)] for (3 = 0.064 and =0.05 (top curve), 0.06, 0.07 and 0.08 (bottom 
curve). For Q > 0.08 we will have more zeros in the interval 1 < r < 10 than 
for Q =0.08. For Q < 0.05 the curves will look like the one for Q = 0.05. 
Then it is easy to see that the integral K(Q) for /(r) = r will be non zero 
for (3 = 0.064 and f2 < 0.06. For n > 1 we will have a more favorable 
situation. In Fig. 4 we plot the integrand of K(Q) for (3=0.064, Q=0.06, 
and f'(r) = df jdr = r (top curve), r 2 /10, and r 3 /100 (bottom curve), i.e., 
for a quadrupolar, octopolar and 16* h -polar external periodic perturbations, 
respectively. In all cases the area under the curve is clearly not null, then 
we will have chaotic motion in these situations. The values of f2= 0.06 and 
0.07 were obtained from fl k ~ 7r/r fc , where corresponds to the maximum 
of the curvature (k = r"/(l + (r') 2 ) 3 / 2 ) of r(r) for the first value, and to 
the inflection point for the second one. From the previous analysis we see 
that the first value of Q, the one based on curvature, is better. Thus, we 
conclude that a perturbation with /(r) ~ r a with a > 1 and Q < Qk will 
always be chaotic. A closer look shows that this criteria is good whenever 
r m ~ r un > 1) otherwise one needs Q <C f^. It is also interesting to compare 
the perturbation frequency Q with uo un = l/r 2 n , that is the angular frequency 
(0) of a particle at the UPO in dimensionless units. For (3 = 0.064 we have 
^«n = 15. So, f2 u un . 

Now we shall study perturbations on the motion of the particle that 
can be modeled by functions f(r) ~ r~ n with n =2, 3,...; n=2 represents 
monopolar contributions, n=3 dipolar, etc. These type of perturbations can 
model forces due to distribution of masses with periodic motions that are 
placed inside the motion. Also, the case /(r) ~ logr, which can be used to 
model the contribution of bodies with a spindly form, will be considered. 

As in the previous case, we shall begin by analyzing a particular situ- 
ation. In Fig. 5 we plot the integrand of K(il) for (3=0.068, f2=0.15, and 
f(r) = l/(3r 2 ) (top curve), l/(9r 3 ), and l/(27r 4 ), (bottom curve). In this 
case we have a rapid oscillation (not shown in Fig. 5) near the point r un 
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(=0.285524) that produces very fine spines of sizes between [-4.0, 4.0] (top 
curve), [-4.5, 4.5], and [-5.5, 5.5] (bottom curve). The areas under the posi- 
tive parts of the curves are clearly greater than the areas under the negative 
ones, which is under very fine spines. So, also in this case K(Q) will be dif- 
ferent from zero and the Melnikov method assures the occurrence of chaos. 
The value £7=0.15 was obtained from £l k ~ It is also instructive to com- 
pare the value of f2 = 0.15 with the UPO frequency. For (3 = 0.068 we find 
u) un = 12. So, again Q <C u un . For powers n>4we have that the relative 
size of the positive part of the area under the curve of the integrand of K (fi) 
begins to decrease and the area of the spines to increase. Hence, the graphic 
method is not appropriated in this case, though one can always numerically 
evaluate K(VL). We will be back to this point later. A similar analysis shows 
that for f(r) ~ logr we also have a strictly positive K(Vt). Therefore, we 
can say that for perturbations characterized by functions f'(r) ~ r~ a with 
1 < a < 4 and O < H^, we will have chaotic orbits. If one further re- 
stricts the range of (3 by the relation r m — r un < 2, one can have chaos with 
perturbations for any value of n > 1 in /(r) ~ r~ n . 

5. Discussion. 

Let us first examine the domain of applicability of our results. The 
method is based in the existence of the homoclinic orbit that restrict the 
parameter (3 in a considerable way (1/16 < (3 < 1/12). With the mass of 
the attraction center k = GM and the areal velocity of the test mass we 
construct our length scale k/h 2 , which is kept arbitrary in our analysis. 

We shall examine as a limit case the value of the parameters for particles 
moving near the minimum of the effective potential v e ff (cf. Fig. 1). In 
principle, these particles will feel to some extent the presence of the saddle 
point (in phase space), the maximum of the potential, that is the responsible 
for the instabilities that give rise to the chaotic motion. In particular, for a 
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circular orbit with radius rp 



Rph 2 /k at the minimum we get, 
JzRqRp 



2(R 



2 4- 
P + 



2 J2R0 



2\2 ' 



(14) 



We have in non dimensional units that r P = (1 + yl — 12p)/2. For the 
characacteristic length of the central body we have r = Roh 2 /k = 2/3/ J 2 . 
Since our potential (0) decays with the distance we need that rp > r . So 
the parameter J 2 that describes the oblateness of the central body plays a 
fundamental role in our analysis . From the limit condition Rq = Rp, we get 
that P = 2 J 2 / (2 + 3 J2) 2 . We recall that J 2 = for a spherical distribution of 
matter, 1/4 for a thin disk, 1/2 for a ring. To be more specific let us consider 
that all the matter is concentrated in a thin disk of radius Rq. We have in 
this case that 0.0625 < (3 < 0.0661. It is instructive to make a table of the 
values of r and rp and e(r ) and e(rp) for different values of (3. 

Table 1. Different parameters for J2 = 1/4 a disk like configuration. 



{3 


r (P) 


rp(P) 


-^( r o) 


-e{r- P ) 


0.0626 


0.708 


.749 


0.593 


1.143 


0.0630 


0.710 


.747 


0.594 


1.142 


0.0640 


0.716 


.740 


0.596 


1.141 


0.0650 


0.721 


.735 


0.599 


1.137 


0.0660 


0.727 


.728 


0.601 


1.136 



Therefore we have some interval of allowed energies and we can have chaos 
without a very fine tuning of the initial conditions. From de values in Table 
1 we conclude that the potential well is rather deep. Furthermore the pa- 
rameters also tells us that we can take for J 2 a value smaller but close to 1/4. 
In this case we can add some structure to the disk: some thickness and/or a 
small central bulge. 

In this work we have presented a graphic method of implementing the 
classic Melnikov method for significant classes of time periodic perturbations 



9 



of a planar motion around an attractive center modeled by the potential (|1|). 
The method, a semi-analytical one, is based on the fact that we were able 
to obtain in a closed form the equation motion of the homoclinic orbit and 
perform a change of variable in the Melnikov function that maps the infinite 
integration interval into a finite one. This allows us to make predictions about 
the zeros of the Melnikov function — therefore about the onset of chaos - 
for significant physically motivated families of periodic perturbations. 

Of course, one can always numerically compute the equivalent function 
K(Q) for any given function f(r) and any value of Q and (3 in the range 
1/16 < (3 < 1/12. Since r(r) is known explicitly, this computation can easily 
be made with a great precision. For the perturbations here considered, /(r) a 
power of r, and a fixed Q > flk, we have that K (Q) takes positive and negative 
values depending on the values of (3. Then for certain specific values of (3 we 
will have K(Q) = and the Melnikov method does not apply in these cases. 
Therefore it does look like that the generic situation is chaotic, i.e., given a 
function f(r) ~ r a and a value of (3, only for a numerable set of frequencies 
Q we should not have chaos. To be more precise, the Melnikov method 
fails to predict chaos only for a numerable set of perturbation frequencies 
fl We note that seldom one can find so easily the parameters range of a 
given problem that will produce a chaotic situation. Therefore, based in our 
previous analysis, we can conjecture that chaos is generic for the class of 
perturbed systems studied in this communication. 

Finally, we want to comment on the fully numerical approaches to the 
Melnikov method. In order to prove numerically the existence of crossings of 



the homoclinic orbits the analyticity of the Hamiltonian is required || [XG 
But numerics can give only indications of the existence of derivatives to a 
finite order. Then the use of the Melnikov method in this case is not in 
conclusive. 

The authors thank CNPq and FAPESP for financial support. 



10 



References 

[1] H. Poincare, Les methodes nouvelles de la mechanique celeste", 
(Gauthier-Villars, Paris 1892). 

[2] H, Bai-lin, Chaos (World Scientific, Singapore, 1982). 

[3] V. I. Arnold Dynamical systems (Encyclopaedia of mathematical sci- 
ences v. Ill, Springer- Verlag, Berlin, 1988). 

[4] V.K. Melnikov, Trans. Moscow Math. Soc. 12, 1 (1963). 

[5] P.J. Holmes, and J.E. Marsden, J. Math. Phys. 23, 669 (1982); ibid, 
Comm. Math. Phys. 82, 523 (1982). 

[6] P. Holmes, Phys. Rep. 193, 137 (1990). 

[7] P.J. Holmes and J.E. Marsden, Indiana U. Math. Jour. 32, 273 (1983). 

[8] O.E. Gerhard,Astr. Asrtrophys. 151, 279 (1985). 

[9] O.E. Gerhard, Mon. Not. R. astr. Soc 222, 287 (1986). 

[10] J. Koiller, J.R.T. de Mello Neto and I. Damiao Soares Phys. Lett. 110A, 
260 (1985). 

[11] L. Bombelli and E. Calzetta. Class. Quantum Grav. 9, 2573 (1992). 

[12] PS. Letelier and W.M. Vieira, Class. Quantum Grav. 14, 1249 (1997). 

[13] R. Moeckel, Comm. Math. Phys. 150 , 415 (1992). 

[14] N. Straumann, General relativity and relativistic astrophysics, 
(Springer- Verlag, Berlin, 1984). 

[15] W.M. Vieira and PS. Letelier, Phys. Rev. Lett. 76, 1409 (1996). 

[16] R.C. Churchill and D.L. Rod, J. Diff. Eq. 37, 23 (1980). 



11 



FIGURE CAPTIONS 

Fig.l. The effective potential is plotted for the dimensionless parameter: 
ft = 0.068 (top curve), 0.072, and 0.078 (bottom curve); the maximum value 
is located at r un = (1 — y/1 — 12/5) /2. 

Fig. 2. A graphic of the positive branch of Eq. (Qj is pictured for values of 
the parameter ft = 0.064 (top curve), 0.065, 0,066, 0,067, and 0.068 (bottom 
curve). 

Fig. 3. A graphic of sin[fh~(r)] for ft = 0.064, and fl =0.05 (top curve), 0.06, 
0.07, and 0.08 (bottom curve). 

Fig.4. Plot of the integrand of K(£l) for (3=0.064, fi=0.06, and f'(r) = 
df /dr = r (top curve), r 2 /10, and r 3 /100 (bottom curve). 

Fig.5. Plot of the integrand of K(Q) for (3=0.068, Q=0.15, and f\r) = 
l/(3r 2 ) (top curve), l/(9r 3 ), and l/((27r 4 ), (bottom curve). 
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